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Abstract 

A computer program ee4f 7 for calculating cross sections of any four fermion final state of e+e~- 
anniliilation at high energy and the corresponding bremsstrahlung reaction that is possible 
in the framework of the Standard Model is presented. As the fermion masses are arbitrary, 
the cross sections for channels that do not contain and/or e~ in the final state can be 
computed without any collinear cut, the on-shell top quark production can be studied and the 
Higgs boson exchange can be incorporated in a consistent way. The program can be used as 
a Monte Carlo generator of unweighted events as well. 



'^Work supported in part by the Polish State Committee for Scientific Research (KBN) under contract No. 
2 P03B 045 23, by the European Community's Human Potential Program under contract HPRN-CT-2000- 
00149 Physics at Colliders and by DFG under Contract SFB/TR 9-03. 

^E-mail: kolodzie@us.edu.pl 

^E-mail: fjeger@ifh.de 



PROGRAM SUMMARY 



Title of program: ee4f7 
Version: 1.0 (February 2004) 
Catalogue identifier: 

Program obtainable from: CPC Program Library or on request by e-mail from the authors 

Licensing provisions: none 
Computers: all 

Operating systems: Unix/Linux 
Programming language used: FORTRAN 90 

CPC Program Library subprograms used: RANLUX 

Cat. Id. Title ref. in CPC 

ACPR RANLUX 79(1994)111 

Memory required to execute with typical data: 4.0 Mb 
No. of bits in a word: 32 

No. of bytes in distributed program, including test data, etc.: 2.016 Mb 
Distribution format: ASCII 

Keywords: e'^e" annihilation, SM, lowest order four fermion reactions, bremsstrahlung, nonzero 
fermion masses 

Nature of physical problem 

Description of all e"^e~ 4 fermions and corresponding bremsstrahlung reactions that are 
possible in the Standard Model (SM) to lowest order and with nonzero fermion masses at 
center of mass energies typical for next generation linear colliders. Such reactions are relevant, 
typically, for W^-pair or intermediate mass Higgs boson production and decay. 

Method of solution 

Matrix elements arc calculated with the hclicity amplitude method. The phase space integra- 
tion is performed numerically utilizing a multi-channel Monte Carlo method. 

Restrictions on complexity of the problem 

No higher order effects are taken into account, except for assuming the fine structure constant 
and the strong coupling at the appropriate scale and partial summation of those one particle 
irreducible loop corrections which are inducing the (fixed) finite widths of unstable particles. 

Typical running time 

The running time depends strongly on a selected channel and desired precision of the result. 
The results of the appended test run have been obtained on a 800 MHz Pentium III processor 
with the use of Absoft FORTRAN 90 compiler in about 9 minutes. In order to obtain a precision 
level below one per mille a few million calls to the integrand are required. This results in less 



than one hour running time for the fastest channel e^e'^ ^pP^i'^^T^r and more than 100 hours 
running time for the slowest channel e^e' —> e^e~e^e~ . However, the typical running time 
for channels that are relevant for the PF-pair or Higgs boson production is several hours. 

LONG WRITE-UP 



1 Introduction 



Precise tests of the Standard Model (SM) that will become possible at future high luminosity 
e+e" linear colliders like TESLA [1^, the NLC [2], or the JLC |31 require adequate high pre- 
cision of theoretical predictions. To meet this requirement is a challenge, as the complexity 
of perturbative calculation increases while going to higher energy and higher luminosity. The 
spectrum of possible reactions involves among others the four fermion channels with possible 
real photon emission. These reactions, which are of particular relevance for W^-pair or inter- 
mediate mass Higgs boson production and decay, receive contributions from several to about 
one thousand Feynman diagrams already at the tree level of the SM. 

In this article, a technical documentation of a numerical program ee4f7 is presented which 
allows for computer simulation of all the four fermion reactions 

e+(pi) + e-{v2) ^ /1(P3) + UV^) + /3(P5) + /4(P6), (1) 

and the corresponding bremsstrahlung reactions 

e+(pi) + e^(p2) ^ /1(P3) + UV^) + /3(P5) + /4(P6) + 7(P7), (2) 

where /j, i = 1, 2, 3, 4 make up a four fermion final state possible in the SM and the particle 
momenta have been indicated in parentheses. The physics contents of the program has been 
discussed in earlier works where the details on the underlying physics, the method of 

calculation and physical results obtained with ee4f 7 have been presented. 

In the approximation of massless fermions, all possible four fermion channels for e^e~ — > 4f 
have been investigated in [7^ (EXCALIBUR) and those for e^e~ 4f7 in j.8j (RacoonWW). ee4f7 
makes possible computations with nonzero fermion masses. Keeping nonzero fermion masses 
will be important in cases when predictions at the 1% accuracy level are required [11-0 • 
Finite masses also provide a natural regularization of distributions which become singular 
in the massless limit. Massive calculations thus provide reliable benchmarks for massless 
calculations with cuts. The latter are much simpler and hence much faster than calculations 
with massive codes. 

ee4f7 may be considered as a building block (the soft plus hard bremsstrahlung part) of a 
complete 0(a) calculation of the processes e^e~ —>■ 4f. Such calculations have been attempted 
in jH] (see also ^]) and have been completed in the double pole approximation for PF-pair 
production [TT] (RacoonWW) and [12] (KORALW/YFSWW) (see also [13^ (EEWW)). The double pole 
approximation incorporates the one-loop corrections for on-shell VT-pair production ^1] and 
the subsequent decay of W^s into fermion-pairs [T3] . 

Several existing programs, some of them being general purpose packages, may be utilized for 
tree-level calculations of e~^e^ 4f, 4f7 with massive fermions. These are GRACE/BASES |16j . 



MADGRAPH/HELAS [HI, CompHEP [IB] (squared matrix element technique), WPHACT NEXT- 
CALIBUR (initial state radiation photons generated via the structure function approach), 
HELAC/PHEGAS [21] (recursive Dyson-Schwinger equation approach), WRAP [22] (ALPHA [23 
algorithm), O'Mega/WH IZARD j21[ and AMEGIC++ [23- Most of the codes work on the basis of 
helicity amplitudes and some use the structure function approach to generate the photons. 
Except for NEXTCALIBUR, which is specialized to e^e~ 4f, all other programs allow to 
generate and evaluate amplitudes for other type of processes. For more details and comparisons 
we refer to the Four Fermion Working Group Report |26j . 

Results obtained with ee4f 7 reproduce nicely, within one standard deviation, the results for 
all channels of (P) and Q in the massless fermion limit of SM that have been surveyed in 
Table 1 of 0. Moreover, results on with nonzero fermion masses of [H] that had been 
obtained with ee4f7 have been reproduced with AMEGIC++ ^7], except for a discrepancy for 
the rather unphysical channels containing on-shell tt-quark-pair, which could be traced back 
to a different implementation of the top quark width in the propagator^. 

The main advantage of the present program is that it allows for relatively fast computation 
cross sections of both and As fermion masses are not neglected, the cross sections 
for channels that do not contain e"^ and/or e~ in the final state can be calculated without 
any kinematical cuts. A number of options have been implemented in the program which 
make possible calculation of the cross sections while switching on and off different subsets of 
the Feynman diagrams. Besides of a possibility of taking into account solely the electroweak 
contributions, or switching off the Higgs boson exchange diagrams, the program allows also 
for a simplified treatment of reactions (Q), (|21) by utilizing two different narrow width approx- 
imations: for W bosons 

e+(pi) + e-(p2) ^ W+{ps,) + W-{p,e) (3) 

^ /1(P3) + /2(P4) + /3(P5) + /4(P6), (4) 

and for a t- quark 

e+(pi) +e-(p2) ^ tips) +t{p456) (5) 

^ t(p3)+&(p4)+/3(P5)+/4(p6). (6) 



2 Scheme of the calculation 



The necessary matrix elements of reactions (Q) and (j21) are calculated with the helicity am- 
plitude method described in |2H| and [Ij and phase space integrations are performed with the 
Monte Carlo (MC) method. 

In order to improve the convergence of the MC integration the most relevant peaks of the 
matrix element squared related to the Breit-Wigner shape of the W, Z, Higgs and top quark 
resonances as well as to the exchange of a massless photon or gluon have to be mapped away. 
As it is not possible to find out a single parametrization of the multi-dimensional phase space 
which would allow to cover the whole resonance structure of the integrand, it is necessary to 

^The misprinted entry for a of e+e^ tbbi at y/s = 500 GeV in Table 1 of [H] should be replaced with 
1.2753(7). 



utilize a multi-channel MC approach [29j. The reader is referred to for more details on the 
multi-channel MC algorithm, on which the MC integration and event generation is based. In 
the next two subsections we collect phase space parametrizations and mappings of integration 
variables that are used in the program. 



2.1 Phase space parametrizations 

The basic phase space parametrizations which are used in the program are described in the 
following. The 4 particle phase space of reaction (P) is parametrized in 2 different ways: 



d^Lips = (27r) ^ d cos 6'34ds34ds56dP5'2 (^34, ml, mfj dPS2 (sse, ml, mfj , (7) 

d'^Lips = (27r)~^ |^ d cos 6'3ds456ds56dPS'2 (^456, S45, mfj dPS2 (sse 



,m^,mQj , 

where s = (pi +P2Y, sij,,, = {pi + pj + ...)^ and dPS'2 (p^, g^, r^) is a two particle (subsystem) 
phase space element defined by 

dPS, (p^ ,^ r^)=SU.p-,-r)'^'^ = j^df!, (9) 

with q being the momentum and fl the solid angle of one of the particles (subsystems) in the 
relative center of mass system, q + f= 0. Momenta ps + pi, ps and the corresponding polar 
angles 6'34, 6^ in Eqs. (jZj) and (jH)) are defined in the centre of mass system (c.m.s.) frame with a 
positive direction of the z axis determined by pi. Using the rotational symmetry with respect 
to the e~^e~ beam axis, the azimuthal angles of ps+pl in Eq. ((7j) and p^ in Eq. (jHJ are set to zero 
which reduces the dimension of the phase space integral from 8 to 7. Parametrization ((7j) covers 
the kinematical situation where the final state fermions make up two 2-particle subsystems 
of approximately equal invariant mass. Similarly, the parametrization (jSI) corresponds to the 
production of a 1- and 3-particle subsystem with the latter decaying again into 1- and 2-particle 
subsystems. 

The 5 particle phase space of reaction Q is parametrized in 3 different ways. Parametrizations 

E 

d^^Lips = (27r)"-^° — dEjd cos 67 dPS'2 (s, S34, s^q) ds34ds56 

X dPS'2 (s34, ml, mfj dPS'2 (-556, ml, mfj , (10) 
E 

d^'^Lips = — dEj dcos6'7dPS'2 (s', S456,m3] ds456ds56 



2 

X dPS2 (s456,S56,m^) dPS'2 (s56,m^, m^) , (11) 

with s' = {pi + P2 — PiY the photon energy E-j and the solid angle ^7, being defined in 
the c.m.s., cover the initial state photon emission. The final state photon emission is covered 
by the following phase space parametrization 

,10 r- _ /o N-io-^^^^("55 "S347, sse) 



d'^'Lips = (27r) 

^ds347ds56d cos 6'347d£'3d£'7dO3d037dPS'2 (s56, ml, mfj , (12) 



8s 

1 , , , 

X 



where the polar angle 6*347 of the momentum P3 + P4 + P7 is defined in the c.m.s.; the energy 
£"3 and spherical angle ^3 of the final state particle that radiates a photon, the photon energy 
£^7 and the azimuthal angle ^37 with respect to the momentum of the radiating particle are 
defined in the frame where P3 + P4 + P7 = 0. Parametrizations (jlUflT^ are also used with 
different permutations of the external particle momenta and we have made use of the rotational 
symmetry with respect to the beam axis in order to reduce the dimension of integration from 
11 to 10. 

The invariant masses Sijk... in Eqs. fl7|)- (jl2|) should, if possible, correspond to the virtuality of 
the propagators of the gauge bosons, Higgs boson and/or top quarks. The peaks related to 
the propagators are then smoothed out by performing suitable mappings pUj . 



2.2 Collinear, soft and quasi real photons 



This subsection is dedicated to the issue of collinear, soft and quasi real photons that is of 
particular relevance for the precision of phase space integration which is performed in the 
program. 

The strong collinear peaking behaviour of the squared matrix element of reaction Q corre- 
sponding to the radiation off the initial state positron 



1/ {Pi-P7 f-ml =-1/ [£7(1-/5 cos 07 



is eliminated by the mapping 



COS07=-^[l-(l+/3e)/r 

Pe 



el ' 



(13) 



(14) 



with re = (1 + /3e)/(l — Pe) and /3e = yl — Am1/s being the velocity of the electron in the 
c.m.s. Similarly, the collinear peaking related to the radiation off the initial state electron 



1/ [(P2 - Pif - 
is dealt with the mapping 



1/[£7(1+/5COS07)] 



(15) 



cos07 = -5-((l-/5e)r:-l) 

Pe 



(16) 



The ~ 1/^7 peaking of the bremsstrahlung photon spectrum of Eqs. (jl3j) and (jl5j) is eliminated 
by the mapping 



nmm 
'7 



E7 = E^''' i^E^-''''/E^ 
where and Ei^^^ are the lower and upper limit of the photon energy £7. 



(17) 



The collinear and soft photon peaking corresponding to radiation off a final state fermion, e.g. 
72(^4), of reaction (j2I) 



1/ iP4+P7Y-ml ~l/(p4-P7)~l/(C3-£3) 



;i8) 



is mapped away with the following mapping 



E, = C,-{Cs-m,)(^ , (19) 

where C3 = y/ 8347/ 2 + {ml - ml)/ (2^5347), with S347 = (^3 + P4 + P?)^- In Eqs- (UnHini) and 
in Eqs. fl24|) - ()28p below, x denotes a random variable uniformly distributed in the interval 
[0,1]. 

In the soft photon limit, the integration over the photon phase space can be performed ana- 
lytically 

I , I 1 /■ d^p7 f Pi P2 , P3 Pa 

l^^^'iprK- = -7^ ^ I ^^i^tt;;: - 3^'^;:::^. + ^^^^ - ^^^^ 

2 5 

d(7o = - I] 9ii9-tjltj, (20) 




where g^i denote the SM couplings of the photon to i-th fermion flavour. The bremsstrahlung 
integrals 7^^, which are defined by 



/S = 7^/ . "w^ , . ■ (21) 



d^P7 Pi ■ Pj 

(27r)3 7|p7|<a; 2^7 {pi ■ pj){pj ■ P7 

for i ^ j may be found in Section 7 of Ref . j^I] . For i = j we have 



where /5i is the velocity of the radiating particle in the c.m.s. and m^ denotes a fictitious mass 
of the photon that has been introduced in order to regularize the infrared divergence. Splitting 
up bremsstrahlung reaction Q into the hard and soft photon parts may serve as a test of 
the quality of integration over the real photon phase space. The inclusive bremsstrahlung 
cross section including both the soft and hard photon parts should be independent of the soft 
photon energy cut u. This independence has been successfully tested for different channels of 
(j21) in The user may repeat such tests on his own by running the program with isoft 

= 1 and ihard = 1 with different values of ecut = u. See Section 3 for explanation of the 
notation. 

Particular attention has also been paid to the channels of (P) and (j21) that include a light 
fermion pair, or an electron and/or a positron in the final state. Matrix elements of such reac- 
tions develop semisingularities at phase space regions where the virtuality of the intermediate 
photon approaches the squared invariant mass of the light fermion pair, or if the photon ex- 
changed in the t-channel becomes quasi-real, as may happen in case of the final state electron 
and/or positron being hardly scattered. 

The decay of the intermediate photon into a light fermion-antifermion pair of invariant mass 
^ff ~ P'jf ~ (P/ is accounted for by the following mapping of the corresponding ~ 1/sjj 

pole in the Feynman propagator 



with s^j^ = {mf + mjY and s™j^ equal to the kinematical upper limit that is computed in 
the program for any specific phase space parametrization. If this virtual photon becomes 



quasi-real it is preferably emitted in the direction parallel to initial beams. The cosine of its 
angle with respect to the beam is then generated with 

cos^ = ?/-(?/ + l)r^', (24) 
in case it is parallel to the e"*" beam, or according to the formula 

cosO = {y -ly-"" -y, (25) 
if it is parallel to the e~ beam. In Eqs. (j24p and (j25|) . 



y 



and r = , (26) 



with Ejj = {Sff + IpjjP)^/^; s' = s for reaction ^ and s' = {pi + P2 — Pi)"^ for reaction (0). 

The pole corresponding to the exchange of a t-channel photon is mapped away by generating 
cosine of the final state electron angle 9^- according to 

cose,- = {y~iy-y, (27) 

and cosine of the final state positron angle 9^- with 

cos^e+ =y-(y + l)/r". (28) 

In Eqs. (EH) and (|2H1), r = (y +!)/(?/- 1) and y = {EiE j - ml) / {\pi\\pf\) , with Ei and Pi {Ef 
and pf) being the energy and momentum of e~ or e"*" in the initial (final) state. 

For the charged current reactions, the integration over cos^e± can be performed to the very 
kinematical limit corresponding to 

t^ = -ml{E,-Eff /{E,Ef). (29) 

However, for some neutral current channels of reactions (^, (0), in particular those including 
another light fermion-antifermion pair in addition to e^e~ in the final state, it may lead to a 
quasi-real pole ~ 1/t in the photon propagator, corresponding to or e~ passing from the 
ininitial to final state with almost no interaction. This singularity can be avoided by imposing 
a cut on the angle of the final state electron or positron with respect to the beam. 



3 Description of the program 

ee4f 7 is a package written in FORTRAN 90. It consists of 93 files including a makefile. They 
are stored in one working directory. The user should specify the physical input parameters in 
module inprms.f and select a number of options in the main program ee4fg.f. 



3.1 Program input 



The default values of the input parameters and options used in the program are those specified 
below. 



3.1.1 Physical parameters 

The initial physical parameters to be specified in module inprms . f are the following. 
Boson masses and widths (GeV): 

• mw=80 . 423 GeV, the W mass, 

• gamw=2. 118 GeV, the W width, 

• mz=91.1876 GeV, the Z mass, 

• gamz=2 . 4952 GeV, the Z width, 

• mh=115.0 GeV, the Higgs mass, 

• gamh=0 GeV, the Higgs width. If gamh = 0, then the Higgs width is calculated according 
to the lowest order prediction of the SM. 

The electroweak (EW) mixing parameter sw2 is then calculated from 

sw2 = l-mw2/mz2 

with 

• mw2 = mw**2 and inz2 = mz**2 in the fixed width scheme, 

• mw2 = mw**2 - i*mw*gainw and mz2 = mz**2 - i*mz*gainz in the complex mass scheme. 
Coupling constants: 

• ralpO = 137 . 03599976 , the inverse of the fine structure constant in the Thomson limit, 

• gmu=l. 16639x10"^ GeV~^, the Fermi coupling constant, 

• alphas=0. 1172, the strong coupling constant at scale mz. 

The inverse of the fine structure constant at scale mw, ralpw, is then calculated from 

ralpw=4 . 44288293815837/ (2*sw2*gmu*mw**2) . 

Fermion masses and widths: 

• me = 0.510998902 MeV, game = MeV, for the electron, 

• mmu = 105 . 658357 MeV, gammu = MeV, for the muon, 

• mtau=l. 77699 GeV, gamtau = GeV, for a tau lepton, 

• mu = 5 MeV, gamu = MeV, for an up quark 

• md=9 MeV, gamd = MeV, for a down quark. 



• mc=1.3 GeV, game = GeV, for a charm quark, 

• ms=150 MeV, gams = MeV, for a strange quark, 

• mt=174.3 GeV, gamt = 1.5 GeV, for a top quark, 

• mb=4 . 4 GeV, gamb = GeV, for a bottom quark. 

No. of colors and conversion constant: 

• ncol=3, the number of colors, 

• convc=0. 389379292x10^^ fb GeV^, a conversion constant. 
3.1.2 Main options 

The following main options should be selected in the main program ee4fg.f. 

• The number of different center of mass (CMS) energies ne 

ne = 1. Recommended if unweighted events are to he generated. 

• The actual values of the CMS energies in the array aecm of size ne: 
aecm=(/500.d0/). 

• The final state of by selecting a value of iproc corresponding to a specific channel 
from the list contained in the file, e.g. 

iproc=l 

corresponds to e"'"e~ ud^~i7^. The bremsstrahlung is then calculated for the same 
channel. 

• Whether or not to calculate the Born cross section, iborn = 1 (yes) / else (no), 
with ncallO calls to the integrand 

iborn = 1 

ncallO = 100000. Recommended No of calls is a few millions. 

• Scan the Born cross section with nscanO calls, iscanO = 1 (yes) /else (no), in order 
to find the dominant kinematical channels, adjust integration weights and find out the 
maximum value of the cross section 

iscanO = 1. This option is strongly recommended. 
nscanO = 2000. A few thousands is recommended. 

• Whether or not to calculate the soft bremsstrahlung cross section, isoft = 1 (yes) / 
else (no), with ncalls calls to the integrand 

isoft = 1 

ncalls = 100000. Recommended No of calls is a few millions. The integration weights 
are the same as those used for the Born cross section. 

• Whether or not to calculate the hard bremsstrahlung cross section, ihard = 1 (yes) 
/ else (no), with ncallh calls to the integrand 

ihard = 1 

ncallh = 100000. Recommended No of calls is a few millions. 



Scan the hard bremsstrahlung cross section with nscanh calls, iscanh = 1 (yes) / else (no) 
iscanh = 1. This option is strongly recommended. 
nscanh = 2000. A few thousands is recommended. 

Generate the unweighted events or not, imc = 1 (yes) /else (no)? 
imc = . 

No standard event record is used. If imc = 1, then the final state particle momenta of 
the accepted unweighted events are printed in file events.dat. 

Include the Higgs boson exchange, ihiggs = l(yes)/else(no)? 
ihiggs = 1 

Calculate the electroweak contributions only, iew = l(yes)/else(electroweaLk + QCD 
Feynman diagrams)? 
iew = 0. 

Whether or not to calculate approximated cross section of (P) itopa = 1 (yes) / else 
(no), with ncalla calls to the integrand, relevant only for top quark pair production 
with an on-shell top quark 
itopa = 

ncalla = 20000. Recommended No of calls is a few hundred thousands. 

Whether or not to calculate approximated cross section of Q iwwa = 1 (yes) / else 
(no), relevant only for W^-pair production 
iwwa = 

This uses an analytic formula and does not require MC integration. 

Choose the scheme: ischeme = Kcomplex mass scheme) /else (fixed width scheme) 
i scheme = 

If ischeme=l, then should alpha_W be complex (iaplw=l) or real (ialpw=0)? 
ialpw=0 . 

3 Auxiliary options 

Choose gauge: iarbg = (unitary gauge) /Karbitrary linear gauge), relevant only 
for the nonradiative top quark production 
iarbg = 

If iarbg=l, then choose gauge parameters > 10^^ corresponds to the unitary gauge): 

ksia=l.dO 

ksiz=l.dl6 

ksiw=l .dl6 

Specify the soft photon energy cut: 
ecut=l .dO 

Impose cuts, icuts > (yes) /else (no)? 
icuts = 1 

If icuts = 1, then specify the kinematical cuts: 



— ctlb = 0.985 - cosine of the charged lepton-beam angle, 

— ctll = 1 - cosine of the angle between charged leptons, 

— ctlq = 1 - cosine of the charged lepton-quark angle, 

— ctgb = . 985 - cosine of the photon-beam angle, 

— ctgl = cos (5 . /180*pi) - cosine of the photon-charged lepton angle, 

— ctgq = cos(5./180*pi) - cosine of the photon-quark angle, 

— ecutg = ecut GeV - minimum hard photon energy, 

— ecutl = 5 GeV - minimum charged lepton energy, 

— ecutq = GeV - minimum quark energy, 

— mqq = 10 GeV - minimum invariant mass of a quark pair, 

— mil = GeV - minimum invariant mass of a charged lepton pair. 

If icuts > 1, then specify the kinematical cuts on a final state electron or positron: 

— theb=5 - minimum angle of a final state electron and/or positron with respect to 
the beam. 

Calculate distributions, idist = 1 (yes)/0(no)? 
idis=0 

If idis = 1, then specify parameters of the distributions: 

nbs = (/ nl , . . . , n8/) - numbers of bins in each distribution. 

The corresponding lower and upper bounds of the distribution arguments 

xmin = (/ xll, . . . , xl8/) - lower bounds, 

xmax = (/ xul, . . . , xu8/) - upper bounds, 

must be specified below in the same file. Constants xli , xui should be of type real (8) 
and ni of type integer. The number of desired distributions and the maximum number 
of bins, nbmajx = ma:x{nl , . . . , n8} should be specified in a module distribs .f . 

The maximum of the fully differential cross section crmax, relevant only if iscanO = 

or iscanh = 0, 

crmax=1000. 



3.2 Routines of ee4fg 

The main program ee4f g, each subroutine, function or module are located in a file named 
the same way as the routine itself, except for subroutines eev, eve and vee, contained in 
gwwv.f, and functions: srr, src and sec, contained in dotprod.f. The package consists of 
the following routines. 

• SUBROUTINE boost - returns a four vector boosted to the Lorentz frame of velocity -v. 

• SUBROUTINE cancuts - checks if generated particle momenta satisfy kinematical cuts. 

• SUBROUTINE carlos - the MC integration routine. 



SUBROUTINE couplsma - returns SM couplings. 
FUNCTION cross - calculates cross sections of and (jSj)- 

FUNCTION crosstopa - calculates approximated cross section of for top quark pair 
production with an on-sliell top quark in the narrow anti-top width approximation. 

The MAIN PROGRAM ee4f g. 

SUBROUTINE cseeww - returns the total unpolarized cross section of on-shell VT-pair 
production. 

SUBROUTINE cutscrit - specifies the class of a process i.e., whether the process is 
purely leptonic, semileptonic, hadronic, etc. 

MODULE cutvalues - contains values of kinematical cuts. 

FUNCTION ddilog - calculates the real part of a dilogarithm for real arguments. 

MODULE distribs - contains parameters of distributions. 

MODULE drivec - contains driving flags and some kinematical variables. 

SUBROUTINE eeee - returns a contraction of a quartic gauge boson coupling with four 
complex four vectors. 

FUNCTION eee - returns a contraction of a triple gauge boson coupling with three complex 
four vectors. 

SUBROUTINE eeff 1 - returns the squared matrix element averaged over initial spins and 
summed over final spins and colors of e~^e~ — >■ //. 

SUBROUTINE eeffffccla - returns the squared matrix element averaged over initial 
spins and summed over final spins and colors of for a charged current processes, 
iproc < 11 or 14 < iproc < 21, and the matrix element squared in the double W 
resonance approximation in arbitrary linear gauge. 

SUBROUTINE eefff f ccl - returns the same as eeffffccla but in the unitary gauge. 

SUBROUTINE eef f f f cc2 - returns the squared matrix element averaged over initial spins 
and summed over final spins and colors of for the mixed charged and neutral cur- 
rent processes, 11 < iproc < 14. These matrix elements can be also obtained with 
eef f f f nc3. 

SUBROUTINE eef f f f cc3 - returns the same as eef f f f cc2 for e^e^ z7gZ/ee"'"e~, iproc = 
21. 

SUBROUTINE eef f f f gidl - returns the squared matrix element averaged over initial spins 
and summed over final spins and colors of the neutral current processes Q containing 
identical leptons, 401 < iproc < 406. 

SUBROUTINE eef f f fgid2 - returns the same as eef f f fgidl for the neutral current pro- 
cesses ((21) containing identical quarks, 406 < iproc < 412. 



• SUBROUTINE eef f f fgidS - returns the same as eef f f fgidl for the neutral current pro- 
cesses (121) containing identical electron neutrinos, iproc = 412. 

• SUBROUTINE eef f f f gncl - returns the squared matrix element averaged over initial spins 
and summed over final spins and colors of the neutral current processes (j2I) containing 
neither virtual gluons nor an external electron, 21 < iproc < 201. 

• SUBROUTINE eef f f fgnc2 - returns the same as eef f f fgncl for the neutral current pro- 
cesses (121) not containing virtual gluons, 201 < iproc < 211. 

• SUBROUTINE eef f f fgnc3 - returns the same as eef ff fgncl for the neutral current pro- 
cesses (12} containing z/gZ/e-pair, 211 < iproc < 301. 

• SUBROUTINE eeffffgncS - returns the same as eeff ffgncl for the mixed charged and 
neutral current processes (j2I), 301 < iproc < 401. 

• FUNCTION eef fffgs - returns the soft bremsstrahlung correction of Eq. ()20|) for any 
channel of (^. 
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• SUBROUTINES eev, eve and vee (contained in gwwv.f) - return contractions of a triple 
gauge boson coupling with two polarization vectors leaving, respectively, the third, sec- 
ond and first Lorentz index uncontracted. 

• SUBROUTINE fefefef - returns matrix elements of the form UifiSFif2SF2TzU2i where 
Fj, i = 1,2,3, and 5*^^, j = 1,2, denote the SM gauge boson-fermion coupling and 
Feynman propagator of a fermion, respectively. 

• SUBROUTINE fefefefh - returns the same as fefefef including contributions of the 
Higgs boson exchange. 

• SUBROUTINE f ef ef - returns matrix elements of the form UifiSFif2U2- 

• SUBROUTINE f ef efh - returns the same as f ef ef including contributions of the Higgs 
boson exchange. 

• SUBROUTINE f ef - returns matrix elements of the form UifiU2- 

• SUBROUTINE f efh - returns the same as f ef including contributions of the Higgs boson 
exchange. 



SUBROUTINE ffvgnew - returns four vectors representing contributions resulting from 
an attachment of a photon to the SM gauge boson-fermion vertex including a possible 
Higgs boson exchange contribution. 

SUBROUTINE f hf - returns matrix elements qhuu^hUU, where Qhuu is the Higgs boson 
coupling to a fermion represented by spinor u and Ah is the Feynman propagator of the 
Higgs boson. 

SUBROUTINE fsfa - returns matrix elements gsuuAs{C)uu, where gsuiU2 is the scalar 
boson coupling to fermions represented by spinors Ui,U2 and A5(^) is the Feynman 
propagator of a scalar boson in an arbitrary linear gauge. 

SUBROUTINE f sf - returns the same as fsfa in the unitary gauge. 

SUBROUTINE fvfa - returns a set of four vectors of the form Ui'y,, 
XM2-Dy'(^), where Ui, U2 are fermion spinors, P± = (1 ± 75)/2, are chirality projec- 
tors, gy^'^ are the gauge boson-fermion couplings of definite chirality and Dy'(^) is the 
propagator of a gauge boson V in arbitrary linear gauge. 

SUBROUTINE f vf - returns the same as fvfa in the unitary gauge. 

MODULE gaugepar - contains driving flags and gauge parameters for calculations in ar- 
bitrary linear gauge. 

MODULE genpvl2 - contains matrix elements representing the photon and Z boson cou- 
pling to initial state fermions. 

MODULE genpvl2g - contains matrix elements representing the photon and Z boson 
coupling to initial state fermions including photon radiation from initial lines. 

MODULE helproc - contains numbers of possible helicity states of external particles. 

SUBROUTINE inipart - calculates content of modules genpvl2 and genpvl2g. 

MODULE inprms - contains initial input parameters. 

FUNCTION Iomega - calculates the soft bremsstrahlung phase space integrals (f2T| and 

(I22D. 

SUBROUTINE kineeff 1 - returns the four momenta, phase space normalization and flux 
factor for a 2 2 process in the CMS. The phase space is parametrized according to 
Eq. ©. 

SUBROUTINE kineef f f f - returns the four momenta, phase space normalization and flux 
factor for a 2 — > 4 process in the CMS. The phase space is parametrized according to 
Eq. ((Zj) or (jHl), dependent on selected options. 

SUBROUTINE kineeffffg - returns the same as kineeff ff for a 2 ^ 5 process in the 
CMS. The phase space is parametrized according to Eq. (fTUjl . |TT|) or (fT^ . dependent 
on selected options. 

SUBROUTINE kineetbf f - returns the same as kineef fff with the phase space parametrized 
according to Eq. (jH)). 



SUBROUTINE kineetbf f g - returns the same as kineef f f f with the phase space parametrized 
according to Eq. (fTT|) . 

SUBROUTINE kinf f - returns the final state four momenta and phase space normahzation 
for a 2 — > 2 process in the CMS. 

SUBROUTINE kinf f g - returns the same as kinf f for a 2 ^ 3 process in the CMS. 
FUNCTION lamsq - the kinematic lambda function, X{^/x, \fz). 
SUBROUTINE mateeee - returns polarized matrix elements of e^e~ — > e~ e^e~e^ . 
SUBROUTINE mateeeeg - returns polarized matrix elements of e+e^ e^e^e^e^7. 
MODULE mathprms - contains necessary arithmetical constants. 
SUBROUTINE matllll - returns polarized matrix elements of e^e~ —>■ llll. 
SUBROUTINE matllllg - returns polarized matrix elements of e^e^ llll'j. 
SUBROUTINE matnnnn - returns polarized matrix elements of e'^e~ z/e^'e^'e^'e- 
SUBROUTINE matnnnng - returns polarized matrix elements of e~^e~ z/gt'ei^e^'eT- 
SUBROUTINE matqqqq - returns polarized matrix elements of e~^e~ qqQQ- 
SUBROUTINE matqqqqg - returns polarized matrix elements of e~^e~ QQQQI- 
MODULE parproc - contains parameters of a specific process. 
MODULE parsm - contains SM couplings and squared masses. 

SUBROUTINE parspec - returns parameters of a specific process ((T)) and calculates 
the top quark and Higgs boson width to the lowest order of SM, the branching ratios for 
the narrow width approximations; initializes weights for the Monte Carlo integration. 

SUBROUTINE psneef f f f - returns a phase space normalization as that of kineef fff for 
a given set of external particle momenta. 

SUBROUTINE psneef fffg - returns a phase space normalization as that of kineef fffg 
for a given set of external particle momenta. 

SUBROUTINE psneetbf f - returns a phase space normalization as that of kineetbf f for 
a given set of external particle momenta. 

SUBROUTINE psneetbf fg - returns a phase space normalization as that of kineetbf fg 
for a given set of external particle momenta. 

SUBROUTINE recpol - returns real polarization vectors of a vector boson in the rectan- 
gular basis. 

FUNCTIONS sec, src and srr (contained in dotprod.f) - return the Minkowski dot 
product of two complex, real and complex, and two real four vectors. 

SUBROUTINE spheric - returns spherical components ps of a four vector p'^, ps= 
( I p I , cos 6, sin 6, cos 0, sin 0) , with 6 and being a polar and azimuthal angles of momen- 
tum p. 



• SUBROUTINE spine - returns the contractions: jPl — p ■ a and + p ■ a, where 'p^ — 

p) is a complex four vector, I is the 2x2 unit matrix and a are the Pauh matrices. 

• SUBROUTINE spinornew ~ returns hehcity spinors in the Weyl representation. 

• SUBROUTINE spinr - returns the same contractions as in spine for a real four vector ]f. 

• MODULE topapprox - contains SM lowest order widths of the top quark in different 
approximations . 

• SUBROUTINE uep - returns a set of spinors of the form u^{g'y^P_ + g'y'' Pj^)Sp^ where Sp 
is a Feynman propagator of an internal fermion and the remaining notation is the same 
as in f vf a. 

• SUBROUTINE usp - returns a set of spinors: u{g^g ^ P_ + g^^^^ P+) S p , where Sp is a Feynman 
propagator of an internal fermion and g^ are complex chirality couplings of the scalar 
boson to the fermion represented by spinor u. 

• SUBROUTINE wwidth - returns the partial width of the 1^-boson averaged over initial 
spins and summed over final spins, and colors for hadronic decay modes. 



3.3 Run output 

A sample of the listing of the test run output is given in an attached file testO. It contains 
a specification of the process, information on the scheme choice, and values of the relevant 
physical parameters in the very beginning. If cuts have been imposed (icut = l) then there 
comes information about specified cuts. Then, if iscan = 1, results of the initial scan and 
resulting weights for the actual run are printed. If calculated weights do not add up exactly 
to 1 the last weight is changed slightly in order to satisfy the normalization condition. This 
change is completely irrelevant numerically, but the message about it is printed. The final 
result for the total cross section is called Integral. Its value is printed together with the 
standard deviation and the actual number of calls used in the calculation. Finally, if imc = 
1 then information on events acceptance efficiency is given, that means a fraction of accepted 
weight 1 events. The same kind of output is printed for the soft (without initial scan) and 
hard bremsstrahlung for each value of the CMS energy. In the very end, all the calculated 
total cross sections together with the corresponding standard deviations are collected in the 
tabular form. 

If the event generation option is switched on, i.e. imc = 1, then the corresponding number 
of unweighted events will be printed out as collections of final state particle four momenta in 
a separate file called events.dat. Whenever a maximum value of the cross section initially 
assigned in ee4fg.f, or found in the result of the initial scan is overfiown a corresponding 
message informing about it is printed in events.dat. The result of integration is still valid, 
however, if the program is run as an event generator, it should be rerun. How to proceed in 
this situation will be described in the next section. 



4 Use of the program 



Up to now the program have been run only on Unix/Linux platforms. In order to run the 
program, the user should select a specific name of a FORTRAN 90/95 compiler, desired options 
and the name of the output file in the makefile. The output file is called test at present. 
The program can be then compiled, linked and run by executing a single command 

make test. 

The results of the test run should reproduce those contained in files testO and events. datO. 

The program can be run as the MC event generator of unweighted events by selecting 
imc = 1 

in ee4fg.f. It is then recommended to run the program for a single CMS energy and to 
perform a scan with a relatively large number of calls, nscanO or nscanh. Attention should be 
paid to possible messages informing about updates of the maximum weight. In this case the 
program should be rerun, however, this time without the initial scan. The initial integration 
weights awO or aw in parspec . f and the maximum weight crmax in ee4f g . f should be updated 
according to the results of the prior scan. At present the efficiency of events acceptance is 
relatively low. However, in view of a relatively fast performance of the program, this should 
not be a serious limitation for the user. 
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